library(flowCore)
library(flowStats)
library(flowViz)
library(ggplot2)
library(lessR)
##install.packages("gridExtra")
library(gridExtra)
##library(multicore)
library(diptest) ## test unimodality

root<-"/space/alpha/dev/FACS/"
##data.dir<-paste(root, "ft_condition_facs", sep="")
data.dir<-paste(root, "ft_condition_2nd_130530", sep="")
setwd(root)


myfplot<-function(wf, data.info, ctl.idx, ncells=10000, xlimi=c(0,800)){
  n.exp<-length(wf)
  par(ps=22, mfrow=c(n.exp,2), mai=c(0.4, 0.4, 0.3, 0.2), oma=c(0,0,0,0))
  ## control 
  xx<-exprs(wf[[ctl.idx]])
  del.idx<-which(xx[,"SSC-H"]<=50)
  xx<-xx[-del.idx,]
  ## select 10000 cells from
  fsc.density<-density(xx[,"FSC-H"])
  fsc.max<-fsc.density$x[which.max(fsc.density$y)]
  ssc.density<-density(xx[,"SSC-H"])
  ssc.max<-ssc.density$x[which.max(ssc.density$y)]
  d<-sqrt((xx[,"FSC-H"]-fsc.max)^2+(xx[,"SSC-H"]-ssc.max)^2)
  o<-order(d)
  sel.idx<-o[1:ncells]
  sel.xx<-xx[sel.idx,]
  cat(fsc.max, ssc.max, "mean:", mean(sel.xx[,"FL1-H"]), "dim:", dim(sel.xx), i, "\n");flush.console()
  plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=19, col="#00000022", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
  points(sel.xx[,"FSC-H"], sel.xx[,"SSC-H"], pch=19, col="#0000ff22", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
  ##plot(density(sel.xx[,"FL1-H"]), col="#000000", ylab="density", xlab="FL1", main=data.info[ctl.idx,23],xlim=c(0,500))
  plot(density(sel.xx[,"FL1-H"]), col="#000000", ylab="density", xlab="FL1", xlim=xlimi, bty="n", main="")
  text(0, mean(density(sel.xx[,"FL1-H"])$y), paste(data.info[1, "effector"], data.info[1, "effector.level"], sep=" "), pos=4)
  for(i in c(1:n.exp)[-ctl.idx]){
    xx<-exprs(wf[[i]])
    del.idx<-which(xx[,"SSC-H"]<=50)
    xx<-xx[-del.idx,]
    ## select 10000 cells from
    fsc.density<-density(xx[,"FSC-H"])
    fsc.max<-fsc.density$x[which.max(fsc.density$y)]
    ssc.density<-density(xx[,"SSC-H"])
    ssc.max<-ssc.density$x[which.max(ssc.density$y)]
    d<-sqrt((xx[,"FSC-H"]-fsc.max)^2+(xx[,"SSC-H"]-ssc.max)^2)
    o<-order(d)
    sel.idx<-o[1:ncells]
    sel.xx<-xx[sel.idx,]
    cat(fsc.max, ssc.max, "mean:", mean(sel.xx[,"FL1-H"]), "dim:", dim(sel.xx), i, "\n");flush.console()
    plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=19, col="#00000022", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
    points(sel.xx[,"FSC-H"], sel.xx[,"SSC-H"], pch=19, col="#0000ff22", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
    plot(density(sel.xx[,"FL1-H"]), col="#000000", ylab="density", xlab="FL1", xlim=xlimi, bty="n", main="")
    text(0, mean(density(xx[,3])$y), paste(data.info[i, "effector"], data.info[i, "effector.level"], sep=" "), pos=4)
  }
}


x<-list()
cnt<-1

file.names<-dir(data.dir, pattern="[0-9]$")
x<-list()
for(i in 1:length(file.names)){
  x[[i]] <- read.FCS(paste(data.dir, "/", file.names[i], sep=""), transformation=FALSE)  
}
wf<-flowSet(x)
##myfplot(wf, sel.data.info, 3)
##wf<-wf[1:14]

ncells<-10000

n.exp<-length(wf)
X11()
nf<-layout(matrix(c(1:30), 5, 6, byrow = FALSE))
layout.show(nf)
##par(ps=22, mfrow=c(5,5), mai=c(0.8, 0.4, 0.3, 0.2), oma=c(0,0,0,0))
par(ps=22, mai=c(0.8, 0.4, 0.3, 0.2), oma=c(0,0,0,0))
## control 

ctl.idx<-c(1,2,2,3,4)
ctl.idx<-c(ctl.idx, c(1,5,5,6,7))
ctl.idx<-c(ctl.idx, c(1,8,8,9,10))
ctl.idx<-c(ctl.idx, c(11,12,13,14,15))
ctl.idx<-c(ctl.idx, c(11,16,17,18,19))
ctl.idx<-c(ctl.idx, c(11,20,21,22,23))

for(i in ctl.idx){
  xx<-exprs(wf[[i]])
  del.idx<-which(xx[,"SSC-H"]<=50 | xx[,"FL1-H"]==0)
  xx<-xx[-del.idx,]
  ## select 10000 cells from
  fsc.density<-density(xx[,"FSC-H"])
  fsc.max<-fsc.density$x[which.max(fsc.density$y)]
  ssc.density<-density(xx[,"SSC-H"])
  ssc.max<-ssc.density$x[which.max(ssc.density$y)]
  d<-sqrt((xx[,"FSC-H"]-fsc.max)^2+(xx[,"SSC-H"]-ssc.max)^2)
  o<-order(d)
  if(nrow(xx)<ncells){
    ncells<-nrow(xx)
  }
  sel.idx<-o[1:ncells]
  sel.xx<-xx[sel.idx,]
  cat(file.names[i], fsc.max, ssc.max, "mean:", mean(sel.xx[,"FL1-H"]), "dim:", dim(sel.xx), i, "\n");flush.console()
  fl1.density<-density(sel.xx[,"FL1-H"])
  fl1.max<-fl1.density$x[which.max(fl1.density$y)]
  ##plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=19, col="#00000022", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
  ##points(sel.xx[,"FSC-H"], sel.xx[,"SSC-H"], pch=19, col="#0000ff22", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
  ##plot(density(sel.xx[,"FL1-H"]), col="#000000", ylab="density", xlim=c(0, max(sel.xx[,"FL1-H"])), xlab="FL1", bty="n", main=file.names[i])
  plot(density(sel.xx[,"FL1-H"]), col="#000000", ylab="density", xlim=c(0,450), xlab="FL1", bty="n", main=file.names[i])
  text(100, 0.003, round(mean(sel.xx[,"FL1-H"]), 2))
}












